Mon. Not. R. Astron. Soc. 000,[T]l5](2002) Printed 30 April 2009 (MN L5TbX style file v2.2) 



IGyr in the Life of the Globular Cluster NGC 6397 
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ABSTRACT 

M4 and NGC 6397 are two very similar galactic globular clusters, which differ mainly in 
their surface brightness profile. M4 has a classic King-like profile, whereas NGC 6397 has a 
more con centrated profi le, which is often interpreted as that of a post-core collapse cluster. 
Heggie & Giersz (|2008j) . however, found that M4 is also a post-core collapse cluster, and 
Giersz & Heggid ( l2009l) concluded that the main reason for the difference between the two 



surface brightness profiles is fluctuations. This conclusion was reached on the basis of Monte 
Carlo models, however, and in the present Letter we verify that similar fluctuations occur in 
A^-body models. The models were initialised by generating initial conditions from the Monte 
Carlo model of NGC6397 at the simulated age of 12Gyr, and one was followed for IGyr. 
The new models help to clarify the nature of the fluctuations, which have the nature of semi- 
regular oscillations with a time scale of order lO^yr. They are influenced by the dynamical 
role which is played by primordial binaries in the evolution of the core. 
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1 INTRODUCTION 

A long-standing problem in the dynamics of galactic globu- 
lar clusters is the observed dichotomy in their surface bright- 
ness profiles. While m ost clusters exhibit a profile similar to 
a classic King profile iKingl 1 19661) . ab out one quarter of the 
, clusters exhibit a more cuspy profile jChemoff & Diorgovskil 
' 1 19891 : IXrager et al.l Il995h . The second class of objects are usu- 
ally described as post-core collapse clusters, as the phenomenon 
of core collapse leads to an object with high central density 
and small core radiu s (Lvnde n-Bell & Wood|[l968l : lLarsorJll97d 

iLvnden-Bell & Eggletoa 198Q , etc. ). By a statistical study of stel- 
lar lum inosity functions, however, lOe Marchi. Paresce, & Pulond 
( l2007h suggested that some clusters with a King-like profile 
might, in fact, be post-core collapse clusters. Independently, 

iHeggie & Giers3 ( l2008t) used a Monte Carlo code for star clus- 
ter evolution to construct a dynamic, evolutionary model of the 
King-like globular cluster M4, and were surprised to find that their 
model was in the post collapse phase of its evolution at the present 
day. They suggested that the reason it exhibited a finite core ra- 
dius is that the core was sustained by heating from a population 
of primordial binary stars. (Many studies had shown that such a 
population, certainly in the case of idealised models with stars 
of equal mass, was sufficient to sustain t he core in this way afte r 
core collapse; see, for examp le, .McMillan. Hut, & Makind ( Il990h : 

iHeggie. Trenti. & H^ ( l2006l) .) 
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The case of the globular cluster NGC 6397 casts serious doubt 
on this interpretation. This cluster has a very similar mass to M4, 
and, if anything, a larger population of primordial boundaries, and 
yet it exhibits a non-King surface brightness profile. If primordial 
binaries sustain the finite core in M4, then NGC 6397 should ex- 
hibit a finite core in the same way. 

This conundrum was considered bv lGiersz & Heggid ( l2009h . 
who constructed a dynamic evolutionary model for NGC 6397, just 
as they had previously done for M4. Not surprisingly, this model 
was also in its post-collapse evolution. Despite having an appro- 
priate primordial binary fraction, however, it exhibited a non-King 
surface brightness profile that was a fair match to that of NGC 
6397. They concluded that primordial binaries and core collapse 
were insufficient to explain the surface brightness profiles of these 
two clusters. By considering and also eliminating other alternatives, 
they concluded that the best explanation for the difference between 
these clusters was one based on fluctuations. They found that, if 
they constructed the surface brightness profile from a single model 
at different times, or from different models started from the same 
initial conditions, but simply with a different seed for the random 
number generator, then the variety of surface brightness profiles 
exhibited differences comparable to the difference between the sur- 
face brightness profiles of M4 and NGC 6397. In effect, their con- 
clusion was the following: that a single cluster could sometimes 
look like M4, and sometimes look like NGC 6397. 

Unfortunately it is not clear that a Monte Carlo code, such 
as the one used for these c luster s, simulates fluctuations correctly. 
ICiersz. Heggie. & HurlevI ( |2008|) did their best to ensure that the 
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Table 1. Data on the initial A'-body model (from a Monte Carlo model at 
12Gyr) 



Number of single stars 


105615 


Number of binaries 


3277 


Total particle number 


112169 


Total mass 


62000 


Mass of binaries 


2400 


Mass of white dwarfs 


27000 


Mass of neutron stars 


3100 


Tidal radius 


22 


Half-mass radius 


3.2 


Core radius 


0.05 



Note: masses are given in A/q, radii in pc. The "core radius" is determined 
as in the Af-body code (Aarseth 2003, pp. 265-266). 

Monte Carlo model behaves similarly to an A'^-body model, at least 
in the range of A'^ where A^-body models are commonly carried 
out, but their investigation was restricted to global quantities, such 
as the evolution of the total mass, and the half-mass radius; they 
did not even consider the question of fluctuations. The nature of 
the fluctuations is also difficult to investigate with a Monte Carlo 
code, which does not follow the orbital motions of the particles. 
These considerations motivate the A'^-body simulations which we 
report in this letter. 

The initial conditions we use are generated from our Monte 
Carlo model of NGC 6397 evolved to the present day, which we 
took to be 12Gyr. We placed the model in a tidal field with a tidal 
radius equal to that of the Monte Carlo model (which used a tidal 
cut-off), but switched off stellar evolution. One of the runs was con- 
tinued for a simulation time of one Gyr. Details and results of the 
runs are given in the following two sections, and the final section 
of the Letter summarises our conclusions. 



2 THE SIMULATIONS 

It is unfortunately impossible to specify the initial conditions with- 
out access to the complete snapshot of the Monte Carlo model of 
NGC 6397 at 12 Gyr, but table 1 gives some summary parameters. 
The Monte Carlo code stores the mass of every star, including bi- 
nary components. The only positional information on a star (or the 
barycentre of a binary) which is held by the Monte Carlo code is 
its radius, and the full position was generated with the assumption 
of spherical symmetry. For the velocity of each star (or barycen- 
tre) the radial and transverse components are available, and the full 
velocity vector was generated by assuming symmetry about the ra- 
dial vector. For a binary, the only internal information (apart from 
the binary masses) is for the semi-major axis and eccentricity. The 
full relative position and velocity were generated assuming a ran- 
dom value of the mean anomaly, and symmetric distributions of the 
orbital plane and the line of apsides. 

The simulation was run with NB0DY6 on a PC equipped with 
a GPU. The CPU was a quad-core Intel Xeon E5410 at 2.33GHz, 
and the G PU a GeForce 980C) GTX . As usual, the code uses A^- 
body units jHeggie & Mathieiilll986h . but the initial virial radius of 
the model (the A'^-body unit of length) was 3.43pc, and its crossing 
time 1.08 Myr (2-\/2 A^-body time units). The half-mass relaxation 
time is of order 700Myr. 

We actually carried out two runs. One was an exploratory (but 
scientifically informative) run with dynamically "inert" binaries, 
i.e. each binary was replaced by a single particle with a mass equal 
to the combined mass of the components. This was run for an equiv- 
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Figure 1. Core and Lagrangian radii (mass fractions 0.1 and 1%) for the 
model with dynamically active binaries (Model A). 
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Figure 2. Core and Lagrangian radius (mass fraction 1%) for the model 
with dynamically inert binaries (Model I). Data on the 0.1% Lagrangian 
radius are not available for this model, and the core radius was stored only 
to 3 decimal places. 

alent of almost 260Myr. The main run used dynamically "active" 
binaries, as described above. We refer to these two runs as "I" 
and "A" respectively. Both simulations proceeded at a rate of about 
IMyr/hr, and the entire run with dynamically active binaries took 
about 1 month for IGyr. 

3 RESULTS AND DISCUSSION 
3.1 Structural evolution 

In the context set out in the Introduction, most of our interest is 
focused on the inner parts of the models, and information on their 
spatial structure is given in Figs [T] and [2] for the models with dy- 
namically active and inert binaries, respectively. The core radius is 
as defined in NB0DY6 ( Aarseth 2003, pp.265-6), and all radii are 
referred to the density centre. 

We note immediately substantial variations in 1 (defined here 
to be the 1% Lagrangian radius). Since the crossing time at this 
radius is of order lO'^yr, it is clear that these are not the kind of 
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fluctuations caused simply by the motion of stars in and out of the 
core, but take place on a much longer time scale. By comparing the 
density within the half-mass radius and 1, and the above estimate of 
the half-mass relaxation time, we may estimate that the relaxation 
time at 1 is under 1 Myr (the unit of time in these figures), and the 
variations of largest amplitude (e.g. the rise between 500 and 700 
Myr in figO take place on a much longer time scale. 

The next obvious observation is that the long-term fluctuations 
have a comparable amplitude for the two runs. But the maximum 
and minimum radii are lower in run I. In particular, the deep min- 
ima in run I are not found in run A. It is natural to attribute these 
facts to the presence of primordial binaries: run I has to make its 
own binaries, which requires higher density than the burning of 
an existing binary. The amplitude of the fluctuations in the inner 
Lagrangian radii is smaller than that of the core radius, but not 
inconsiderable, resulting in variations in the mean density within 
the 1% Lagrangian radius by a factor of around 3:1. 

Though we have not shown data on the outer Lagrangian radii, 
we find that the variations we observe are largely confined to the 
inner few percent of the mass. In fact the 10% Lagrangian radius 
fluctuates with a relative amplitude of order 7%, and the oscillations 
are approximately 180° out of phase with those of the core. The 
relative amplitude of oscillations drops to less than 1% at the half- 
mass radius. 

The initial decrease in both runs is an interesting feature. One 
might suspect that it is due to some property of the Monte Carlo 
code, which perhaps maintains thermal equilibrium in a structure 
which would not be in thermal equilibrium in an A^-body model. 
Another possibility is a response of the system to the cessation of 
mass loss by stellar evolution. And yet, in the overall range of data 
in these figures, the structure at the start is not unusual. Further- 
more, it must be remembered that this model was selected because, 
at 12 Gyr, it yielded a surface brightness profile resembling that of 
NGC6397, and this happens only intermittently, even in the Monte 
Carlo model. 

Fig[3] when considered in conjunction with Figs[T]and[2l al- 
lows a direct qualitative comparison between the Monte Carlo and 
A'^-body models in respect of the fluctuations in the inner radii, 
though the definition of the core radius adopted in the Monte Carlo 
model is different (see caption). This figure shows oscillations of 
a similar amplitude and range of time scales to those exhibited by 
Run A (FiglB. 

In an attempt to make this statement more quantitative we 
show in Fig|4] the autocorrelation of the core radius for Run A 
and the Monte Carlo model, though it should be borne in mind 
that these correspond to different time intervals. In both cases we 
have restricted the offset r in the autocorrelation (defined to be 
{x{t)x(t + r)) for a zero-mean normalised signal x{t)) to be less 
than half the duration of the measurements. The result for the A*'- 
body model is striking, and confirms the visual impression from 
FiglUof fairly regular oscillations with a period of about 400Myr. 
There is a faint suggestion of similar structure in the autocorrelation 
function of the Monte Carlo model. 

Another structural quantity of interest is, of course, the to- 
tal mass, and the A'^-body run provides a check on this aspect of 
the Monte Carlo model. In the A'^-body model the rate of mass 
loss increases for the first few tens of Myr, presumably because 
escapers take some time to reach the boundary (at 2 tidal radii) 
where escapers are removed. For the period after 200 to 400Myr, 
however, the rate of mass loss is virtually constant, and yields 
d{\nM)/dt = -1.35X 10"*±2x lO"'^, wheretheunitof timeis 
IMyr. For the Monte Carlo model discussed above, we do not have 
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Figure 3. Core and Lagrangian radius (mass fraction 1%) for the Monte 
Carlo model (which has dynamically active binaries) between 10 and 12 
Gyr. Data on the 0.1% Lagrangian radius are not available for this model, 
and data were stored with lower frequency. The core radius Tc is defined 
t'y ''c = 3(''''^)/(4vrGpo), where the mass-weighted mean square three- 
dimensional velocity {v'^) and the central density po are calculated for the 
innermost 20 stars. 
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Figure 4. Autocorrelation of the core radius in Run A and in the Monte 
Carlo model. For the latter, which was output at irregular intervals, the data 
was interpolated. The abscissa r is defined in the text. 



data extending beyond 12Gyr. For another model differing only in 
the initial seed, however, in the same period the corresponding re- 
sult is d(ln M) /dt = -1.546 x 10"* ± 5 x 10"'^. This difference 
of 13% should be corrected for the fact that there is no mass loss 
through stellar evolution in the A'-body model. In the Monte Carlo 
model this contributes only about 3% of the total, but the direct 
loss of this mass also induces further loss of stars (by tidal over- 
flow) because it make the potential well more shallow. It is difficult 
to quantify this induced mass loss, but it appears that the discrep- 
ancy in the total rate of mass loss between the two models is less 
than 10%. 
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Surface density profiles 
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Figure 5. Surface densities of three groups of stars in Run A at two times. 
Only the central 30 arcsec (approximately) is shown. For neutron stars the 
surface density is reduced by a factor 10, to avoid confusion with the suiface 
density of white dwarfs. 



3.2 Surface brightness profiles 

Since so much of the mass at small radii is in the form of degen- 
erate remnants, the influence of the fluctuations shown in Fig[T]on 
the surface brightness profile is not obvious. Insufficient data were 
collected from Run I, but Fig[5] provides a useful comparison of 
two snapshots from Run A. These are taken at times of 200 and 
340 Myr, which correspond, respectively, to low and high values of 
the radii plotted in Fig[T] (though not the extreme values). As ex- 
pected, the central surface densities are lower at the time when the 
Lagrangian and core radii are larger. 

Let us consider these profiles in more detail. Neutron stars 
and white dwarfs are the dominant components at the centre, each 
contributing about half of the projected number density. (Note that 
the surface density of neutron stars in Fig[5]is divided by 10, to 
avoid overlap with the white dwarfs.) Beyond the core, however, 
the projected density of neutron stars decreases more strongly, as 
is expected from their greater mass. Indeed they are actually the 
dominant central component in terms of spatial mass density. The 
component which contributes least to the projected central density, 
by a factor of order 10, are the non-degenerate stars, though they 
begin to dominate at radii beyond those shown in the figure. 

All three components show a lower surface density and a 
larger core at the later epoch (340 Myr) compared to the earlier 
epoch (200 Myr). It is interesting to attempt to quantify this in terms 
as close as possible to observational procedures, by considering the 
radius at which the projected density of the non-degenerate com- 
ponent drops below its central value by a factor of two. Despite the 
spatial fluctuations, it is clear that this radius is about 15" at the 
later epoch, and about 5" at the earlier one. Actually, it would also 
be reasonable to describe the latter profile as a shallow cusp. In- 
deed, in view of the very small number of stars within 1 arcsec, this 
is a more robust description. 

At radii beyond about 4 arcsec the changes in the surface den- 
sity are much smaller but go in the opposite direction, at least for 
the degenerate components. At the epoch when the core and La- 
grangian radii are larger, the number density at these radii is higher, 
the total mass of each component being conserved (except for es- 
cape). 

This discussion still does not necessarily correspond to an 
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Figure 6. Number and internal energy of b inarie s in run I. The units of 
energy are A^-body units >Heggie & Mathiel3ll986l) . 



observational description, which is usually expressed in terms of 
surface brightness, i.e. weighted by luminosity. In particular, the 
most luminous stars are centrally concentrated relative to most non- 
degenerate stars. Because stellar evolution is absent from the A*'- 
body simulation, however, it is not an entirely appropriate model 
for considering how fluctuations affect the surface brightness. T he 
discussion of our Monte Carlo model in lOiersz & Heggi3 ( l2009h is 
more complete in this respect. 



3.3 Binaries 

Run A contains over 3000 binaries, and it is best to begin with 
run I, where the effects of individual binaries are clearer (Fig[6l(. 
By comparing with Fig[2l we can surmise that one or two bina- 
ries were responsible for terminating the core collapse at about 
70Myr and initiating the subsequent expansion. For a period be- 
ginning at I25Myr, however, no binaries were present, suggesting 
that the continuing modest exp ansion was powered gravothermally 
dSugimoto & Bettwiesej|T983h . A second phase of core collapse 
and binary formation appears to have started at about 225Myr, but 
is incomplete by the end of the run. 

Fig|7] shows similar data for run A. The numbers are not in- 
cluded, but the binaries are separated into two groups. Those la- 
belled "new" exclude binaries which have evolved from primordial 
binaries by exchange. If, however, two binaries collide, resulting 
in a hierarchical triple system, the outer motion of the hierarchy 
is regarded as a "new" binary. In any event we see that both types 
of binary have an active role to play. The correlation with the evo- 
lution of the radii (FiglTJ is less clear than for run I, though the 
initial contraction seems to be associated with a period of relatively 
sluggish binary activity. 

While neither run shows the kind of gravothermal oscillations 
which are so evident in simulations of systems with equal masses, 
it should not be surprising to find evidence of gravothermal effects 
in a multi-mass system of the kind which we are studying (see the 
discussion in Giersz & Heggie (2009)). 

The rate of change of the binary fraction in run A is +2.05 x 
10~^ ± 2.5 X 10~*, where the unit of time is IMyr. For the Monte 
Carlo model the corresponding value is+1.56xl0~^±L5xl0~*, 
and these results are pretty consistent within the errors. 
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Figure 7. Internal energy of primordial and new binaries in run A. Though 
the origins differ (because of the enormous energy of the primordial binary 
population) the scales on the vertical axes are equal. The meaning of "new" 
and "primordial" is discussed in the text. 



4 CONCLUSIONS 

We have carried out an AT-body simulation of the globular clus- 
ter NGC6397, in order to study the evolution of its central struc- 
ture over a period of IGyr starting at the present day. The sim- 
ulation was initialised u sing the results of a Monte Carlo model 
jGiersz & Hegri3l2009l) which approximately fits the present-day 
profiles of surface brightness and velocity dispersion, and the mass 
function at two radii. The main limitation of the A'^-body model is 
that there is no stellar evolution. This apart, it is in many respects 
the most reaUstic A^'-body simulation of a specific globular cluster 
of which we are aware. 

The model provides a new check on the reliability of the 
Monte Carlo code, and suggests that the evolution of the binary 
fraction is satisfactory, while the rate of escape of mass appears to 
be in agreement to better than 10%. 

The results show that the population of primordial binaries 
in the cluster (about 3%) suppresses the deepest collapses of the 
core, which nevertheless still exhibits substantial fluctuations on a 
time scale of many core relaxation times. Their amplitude is suffi- 
cient to change the mean density of the innermost 1% of the mass 
by a factor of order 3:1. These changes are reflected in variations 
in the core radius, whether measured in terms of the spatial den- 
sity, as in an A'^-body model, or by the radius at which the surface 
density decreases to half its central value. These changes manifest 
themselves in all components that we have studied: neutron stars 
(the dominant component in the central density), white dwarfs, and 
non-degenerate stars, which are the least dominant. 

In view of our results, it is interesting to think of globular clus- 
ters as variable stellar systems, in much the same way that many 
stars are variable. The time scales and mechanisms are vastly differ- 
ent, and the variations in globular clusters are confined to the vicin- 
ity of the core. But it is another demonstration of the deep physical 
resemblance between these two types of thermal, self-gravitating 
objects. 
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